Analyzing and summarizing the effects of time on exercise in the murine adipose tissue

---

Packages

rm(list = ls())

easypackages::libraries("edgeR",

                        "data.table",
                        "dplyr",
                        "org.Hs.eg.db",
                        "org.Mm.eg.db",
                        "openxlsx",
                        
                        "ggplot2",
                        "patchwork",
                        "ggsci",
                        
                        "reshape2",
                        
                        "UpSetR",
                        "clusterProfiler",
                        "ReactomePA",
                        "clipr",
                        "beepr") |>
  suppressPackageStartupMessages()

## Warning: package 'edgeR' was built under R version 4.1.1

## Warning: package 'S4Vectors' was built under R version 4.1.1

## Warning: package 'clusterProfiler' was built under R version 4.1.1

## Warning: package 'clipr' was built under R version 4.1.2

## All packages loaded successfully

upset <- function(labjust = -40, ...) {
  old <- theme_set(theme_grey() + theme(axis.title.y.left = element_text(vjust = labjust)))
  print(UpSetR::upset(...))
  # UpSetR::upset(...)
  theme_set(old)
}

theme_set(theme_bw() +
            theme(axis.title = element_text(size = 8, color = "black"),
                  axis.text = element_text(size = 5, color = "black"),
                  plot.title = element_text(size = 10, color = "black")))

Introduction

This script summarizes and shows the code for all figures generated for the project “Time-of-day influences post-exercise metabolism in mouse adipose tissue”. RNA sequencing has been performed by Lucile Dollet and Logan Pendergast.

Analysis

Load data

load("data/data.all.Rdata")

dataDGE <- DGEList(counts = data.raw, 
                    remove.zeros = T,
                    samples = dataGroups, 
                    genes = rownames(data.raw))

## Removing 13508 rows with all zero counts

dataDGE$genes$Symbol <- mapIds(org.Mm.eg.db, 
                          dataDGE$genes$genes,
                          keytype="ENSEMBL", 
                          column="SYMBOL")

## 'select()' returned 1:many mapping between keys and columns

dataDGE$genes$entrez <- mapIds(org.Mm.eg.db, 
                          dataDGE$genes$genes,
                          keytype="ENSEMBL", 
                          column="ENTREZID")

## 'select()' returned 1:many mapping between keys and columns

dataDGE <- dataDGE[!is.na(dataDGE$genes$Symbol),]
rm(data.raw)

Quality control

Distribution at various filters

Filter by group to capture things that are high in one group and low in another (instead of average across all samples)

design <- model.matrix(~0 + super_group, data = dataGroups)
colnames(design) <- gsub(colnames(design), pattern = "super_group", replacement = "")

g <- lapply(c(0,5:15), \(i){
  keep <- filterByExpr.DGEList(dataDGE, design = design, min.count=i)
  temp <- dataDGE[keep,,]
  temp <- calcNormFactors(temp)
  temp <- cpm(temp,log = T, prior.count = 0.25) |> as.data.table()
  temp <- data.table::melt(temp,stringsAsFactors = F) |> suppressWarnings()#melt warning
  temp <- cbind(temp, dataGroups[match(temp$variable, dataGroups$basespaceID), 
                                 c("basespaceID", "exercise_group", "biopsy_time")])

  temp <- ggplot(temp, aes(x = value, group = basespaceID)) +
    geom_density(adjust=1, show.legend=F) +
    facet_grid(exercise_group ~ biopsy_time) +
    ggtitle(paste(i," transcript filter density",sep = ""))
  
  temp
})
names(g) <- paste0("filt_",c(0,5:15))
g$filt_0 + g$filt_4 + g$filt_10 + g$filt_15

Filter by 7 looks as good as filter by 15. Going with 10.

Filtering by expresion and group

design <- model.matrix(~0 + super_group, data = dataGroups)
keep <- filterByExpr.DGEList(dataDGE, design = design, min.count = 10) #10 is default

dataDGE <- dataDGE[keep,,keep.lib.size=F]
dataDGE <- calcNormFactors(dataDGE)
dataDGE <- estimateDisp(dataDGE, design)

Export for GEO submission

out <- cbind(dataDGE$genes,cpm(dataDGE))
colnames(out)[1] <- "Ensembl"
rownames(out) <- NULL
fwrite(out, "GEO/GEO upload/processed.data.csv")

MDS plot with all variables (post filter).

Ploting using plotMDS.DGElist (code not shown)

Dim 1 & 2

par(mfrow=c(2,2))
plotMDS.DGEList(x = dataDGE, col = as.numeric(as.factor(dataDGE$samples$exercise_group)), labels = dataDGE$samples$samples)
title("exercise_group")
plotMDS.DGEList(x = dataDGE, col = as.numeric(as.factor(dataDGE$samples$SeqPool)), labels = dataDGE$samples$samples)
title("SeqPool")
plotMDS.DGEList(x = dataDGE, col = as.numeric(as.factor(dataDGE$samples$biopsy_time)), labels = dataDGE$samples$samples)
title("biopsy_time")
plotMDS.DGEList(x = dataDGE, col = as.numeric(as.factor(dataDGE$samples$super_group)), labels = dataDGE$samples$samples)
title("super_group")

The samples that are kinda sticking out in the mds relativly low read counts, but all are above 24.9 million.

Differential expression

dataGroups$super_group <- factor(dataGroups$super_group, levels = c("SED_m", "SED_e", "RUN_m", "RUN_e"))
dataDGE$samples$super_group <-  factor(dataDGE$samples$super_group, levels = c("SED_m", "SED_e", "RUN_m", "RUN_e"))

design <- model.matrix(~0 + super_group, data = dataDGE$samples)
colnames(design) <- gsub(colnames(design), pattern = "super_group", replacement = "")

fit <- glmQLFit(dataDGE, design)

colnames(design)

## [1] "SED_m" "SED_e" "RUN_m" "RUN_e"

contrasts <- makeContrasts(
  run_evening_vs_moring = RUN_e - RUN_m,
  sed_evening_vs_moring = SED_e - SED_m,
  morning_exercise = RUN_m - SED_m,
  evening_exercise = RUN_e - SED_e,
  levels=design)

res <- NULL

for(i in colnames(contrasts)){
  qlf <- glmQLFTest(fit, contrast=contrasts[,i])
  res[[i]] <- topTags(qlf, n = Inf, sort.by = "none")$table |>
    as.data.table()
}

resSig <- lapply(res, function(x){
  x <- x[x$FDR<0.05,]
  x
})

counts <- cpm(dataDGE, log = T, prior.count = 0.25)
rownames(counts) <- dataDGE$genes$Symbol

Checking p-value distributions

par(mfrow = c(2,2))
lapply(names(res), function(x) print(hist(res[[x]]$PValue, main = x)))

Summarizing results

Upset

dataUpset <- lapply(X = resSig, function(x){
  x$genes
})

names(dataUpset) <- c("Time of day effect exercise",
                      "Time of day effect sedentary",
                      "Rest phase exercise",
                      "Active phase exercise")
png("figures/upset.png", height = 11.69/2, width = 8.27, units = "in", res = 300)
upset(labjust = -30,
      fromList(dataUpset), 
      text.scale = 1.5, 
      point.size = 5,
      set_size.angles = 0)
dev.off()

## quartz_off_screen 
##                 2

upset(labjust = -30,
      fromList(dataUpset), 
      text.scale = 1.5, 
      point.size = 5,
      set_size.angles = 0)

LogFC

Comparing logFC contrasts

Circadian exercise effect

dataPlot <- data.table(SYMBOL = res$run_evening_vs_moring$Symbol,
                       Sed = res$sed_evening_vs_moring[,c("logFC", "FDR")],
                       Run = res$run_evening_vs_moring[,c("logFC", "FDR")])

sig <- function(x, cont1, cont2){
  cont1test <- as.numeric(x[paste0(cont1,".FDR")]) < 0.05
  cont2test <- as.numeric(x[paste0(cont2,".FDR")]) < 0.05
  if(cont1test & cont2test)
    return("Both")
  else if(cont1test)
    return(cont1)
  else if(cont2test)
    return(cont2)
  else
    return("None")
}

dataPlot$sig <- apply(dataPlot, 1, sig, cont1 = "Sed", cont2 = "Run") |> as.factor()
levels(dataPlot$sig)[3:4] <- c("Exercised", "Sedentary")

dataPlot$highlight <- "no"

dataPlot[xor(abs(Sed.logFC)>2.5,
             abs(Run.logFC)>2.5) & 
           sig %in% c("Exercised", "Sedentary", "Both"),
         highlight :="yes"]

dataPlot[(abs(Sed.logFC)>2.5 |
           abs(Run.logFC)>2.5) & 
           sig %in% c("Exercised", "Sedentary", "Both"),
         highlight :="yes"]

#removing some overploting running genes
dataPlot[sig == "Exercised" & Run.logFC>-3 & Run.logFC < 0,
         highlight :="no"]

ggplot() +
  geom_vline(xintercept = 0, color = "gray") +
  geom_hline(yintercept = 0, color = "gray") +
  geom_abline(slope = 1, color = "gray", lty = "dashed") +
  geom_point(data = dataPlot[sig == "None"],
             aes(x = Run.logFC, y = Sed.logFC),
             color = alpha("gray", .2),
             size = 1) +
  geom_point(data = dataPlot[sig != "None"],
             aes(x = Run.logFC, y = Sed.logFC),
             color = "black",
             size = 1.1) +
  geom_point(data = dataPlot[sig != "None"], 
             aes(x = Run.logFC, y = Sed.logFC, color = sig),
             size = 1) +
  geom_text(data = dataPlot[highlight == "yes"], 
                  aes(x = Run.logFC, y = Sed.logFC, label = SYMBOL), 
                  nudge_y = 0.3,
                  size = 1.5) +
  scale_color_igv(name = "Significance") +
  xlab("Time of day effect in exercised mice (logFC)") +
  ylab("Time of day effect in sedentary mice (logFC)") +
  theme(legend.text = element_text(size = 5),
        legend.title = element_text(size = 7))

ggsave("figures/Time of day effect.png", height = 297/4, width = 210/2, units = "mm")

dataPlot <- data.table(SYMBOL = res$run_evening_vs_moring$Symbol,
                       "Rest phase" = res$morning_exercise[,c("logFC", "FDR")],
                       "Active phase" = res$evening_exercise[,c("logFC", "FDR")])
dataPlot$sig <- apply(dataPlot, 1, sig, cont1 = "Rest phase", cont2 = "Active phase") |> as.factor()

dataPlot$highlight <- "no"

dataPlot[xor(abs(`Rest phase.logFC`)>2.5,
             abs(`Active phase.logFC`)>2.5) & sig %in% c("Rest phase", "Active phase"),
         highlight :="yes"]
dataPlot[highlight == "yes" & nchar(SYMBOL)>=7, highlight := "no"]

ggplot() +
  geom_vline(xintercept = 0, color = "gray") +
  geom_hline(yintercept = 0, color = "gray") +
  geom_abline(slope = 1, color = "gray", lty = "dashed") +
  geom_point(data = dataPlot[sig == "None"],
             aes(x = `Active phase.logFC`,y = `Rest phase.logFC`),
             color = alpha("gray", .2),
             size = 1) +
  geom_point(data = dataPlot[sig != "None"],
             aes(x = `Active phase.logFC`,y = `Rest phase.logFC`),
             color = "black",
             size = 1.1) +
  geom_point(data = dataPlot[sig != "None"], 
             aes(x = `Active phase.logFC`,y = `Rest phase.logFC`, color = sig),
             size = 1) +
  geom_text(data = dataPlot[highlight == "yes"], 
            aes(x = `Active phase.logFC`,y = `Rest phase.logFC`, label = SYMBOL), 
            nudge_y = 0.2,
            size = 1.5) +
  scale_color_manual(name = "Significance", values = "#e7298a") +
  xlab("Active phase exercise (logFC)") +
  ylab("Rest phase exercise (logFC)") +
  scale_x_continuous(limits = c(-5,5), breaks = seq(-5,5, by = 2.5)) +
  scale_y_continuous(limits = c(-5,5), breaks = seq(-5,5, by = 2.5)) + #warning for removing 3 nonsig
  theme(legend.text = element_text(size = 5),
        legend.title = element_text(size = 7))

## Warning: Removed 13 rows containing missing values (geom_point).

ggsave("figures/Exercise effect.png", height = 297/4, width = 210/2, units = "mm")

## Warning: Removed 13 rows containing missing values (geom_point).

Individual genes

Why are some with large effect size not sig? High variance obviously.

targets <- c("Acnat2", "Fos", "Nr1d1", "Nr4a3", "Tchh", "Ighd")

dataPlot <- reshape2::melt(counts[targets,])
dataPlot <- merge(dataPlot, data.frame(dataDGE$samples, Var2 = rownames(dataDGE$samples)), by = "Var2")
dataPlot$cpm <- 2^dataPlot$value
dataPlot$biopsy_time <- as.factor(dataPlot$biopsy_time)
levels(dataPlot$biopsy_time) <- c("ZT15", "ZT3")
dataPlot$biopsy_time <- relevel(dataPlot$biopsy_time, "ZT3")
dataPlot$exercise_group <- as.factor(dataPlot$exercise_group)
levels(dataPlot$exercise_group) <- c("Exercised", "Sedentary")

plotGenes <- function(){
  ggplot(dataPlot, 
         aes(x = super_group, y = cpm, color = biopsy_time, shape = exercise_group)) +
    ggbeeswarm::geom_quasirandom(width = 0.25, size = 1.5, color = "black") +
    ggbeeswarm::geom_quasirandom(width = 0.25, size = 1.25) +
    stat_summary(fun = "mean", geom = "crossbar", inherit.aes = T, color = "black", size = 0.1) +
    scale_y_log10() +
    scale_x_discrete(labels = c("Sedentary\nZT3", "Sedentary\nZT15","Exercised\nZT3", "Exercised\nZT15")) +
    xlab("") + ylab("Expression (cpm)") +
    facet_wrap(~Var1, scales = "free_y") +
    scale_color_manual(name = "Time of day", values = c("#FE711C","#0707D6")) +
    theme(strip.background = element_rect(fill = "gray95"),
          legend.position = "bottom",
          legend.text = element_text(size = 6),
          legend.title = element_text(size = 6),
          legend.margin = margin(-20,0,0,0, unit = "pt"))
}

plotGenes()

ggsave("figures/Selected sig genes.png", height = 297/4, width = 210*0.66, units = "mm")

Visualizing the high logFC that are not significant at the left most quadrant of the time of day effect

targets <- intersect(res$run_evening_vs_moring[logFC<(-5.5)& FDR>0.05, Symbol],
                     res$sed_evening_vs_moring[logFC>(4) & FDR>0.05, Symbol])

dataPlot <- reshape2::melt(counts[targets,])
dataPlot <- merge(dataPlot, data.frame(dataDGE$samples, Var2 = rownames(dataDGE$samples)), by = "Var2")
dataPlot$cpm <- 2^dataPlot$value
dataPlot$biopsy_time <- as.factor(dataPlot$biopsy_time)

levels(dataPlot$biopsy_time) <- c("ZT15", "ZT3")
dataPlot$biopsy_time <- relevel(dataPlot$biopsy_time, "ZT3")

dataPlot$exercise_group <- as.factor(dataPlot$exercise_group)
levels(dataPlot$exercise_group) <- c("Exercised", "Sedentary")

plotGenes()

ggsave("figures/Selected non sig genes time of day.png", height = 297/4, width = 210*0.66, units = "mm")

Visualizing the high logFC that are not significant at the left most quadrant of the time of exercise effect

targets <- intersect(res$evening_exercise[logFC<(-4.5)& FDR>0.05, Symbol],
                     res$morning_exercise[logFC>(2) & FDR>0.05, Symbol])[1:6]

dataPlot <- reshape2::melt(counts[targets,])
dataPlot <- merge(dataPlot, data.frame(dataDGE$samples, Var2 = rownames(dataDGE$samples)), by = "Var2")
dataPlot$cpm <- 2^dataPlot$value
dataPlot$biopsy_time <- as.factor(dataPlot$biopsy_time)
levels(dataPlot$biopsy_time) <- c("ZT15", "ZT3")
dataPlot$biopsy_time <- relevel(dataPlot$biopsy_time, "ZT3")

dataPlot$exercise_group <- as.factor(dataPlot$exercise_group)
levels(dataPlot$exercise_group) <- c("Exercised", "Sedentary")

plotGenes()

ggsave("figures/Selected non sig genes exercise.png", height = 297/4, width = 210*0.66, units = "mm")

Enrichment

goMF <- compareCluster(gene = lapply(resSig[c("morning_exercise",
                                              "evening_exercise")], \(x) x$entrez),
                       universe = res$run_evening_vs_moring[, entrez], 
                       OrgDb = org.Mm.eg.db,
                       ont = "MF",
                       fun = "enrichGO",
                       readable = T)
goMF

## #
## # Result of Comparing 2 gene clusters 
## #
## #.. @fun      enrichGO 
## #.. @geneClusters    List of 2
##  $ morning_exercise: chr(0) 
##  $ evening_exercise: chr [1:121] "20671" "20353" "110895" "67876" ...
## #...Result   'data.frame':   7 obs. of  10 variables:
##  $ Cluster    : Factor w/ 2 levels "morning_exercise",..: 2 2 2 2 2 2 2
##  $ ID         : chr  "GO:0001228" "GO:0001216" "GO:0001968" "GO:0003707" ...
##  $ Description: chr  "DNA-binding transcription activator activity, RNA polymerase II-specific" "DNA-binding transcription activator activity" "fibronectin binding" "steroid hormone receptor activity" ...
##  $ GeneRatio  : chr  "13/107" "13/107" "4/107" "4/107" ...
##  $ BgRatio    : chr  "356/14341" "360/14341" "31/14341" "43/14341" ...
##  $ pvalue     : num  2.41e-06 2.73e-06 7.89e-05 2.89e-04 2.89e-04 ...
##  $ p.adjust   : num  0.000378 0.000378 0.007287 0.013341 0.013341 ...
##  $ qvalue     : num  0.00036 0.00036 0.00695 0.01272 0.01272 ...
##  $ geneID     : chr  "Sox17/Nr4a2/Nr4a3/Fosl2/Zfp606/Zfp930/Dlx3/Sox4/Irf4/Maff/Nr4a1/Litaf/Zfp948" "Sox17/Nr4a2/Nr4a3/Fosl2/Zfp606/Zfp930/Dlx3/Sox4/Irf4/Maff/Nr4a1/Litaf/Zfp948" "Thbs1/Sdc4/Epha1/Vegfa" "Nr4a2/Nr4a3/Vdr/Nr4a1" ...
##  $ Count      : int  13 13 4 4 4 4 3
## #.. number of enriched terms found for each gene cluster:
## #..   morning_exercise: 0 
## #..   evening_exercise: 7 
## #
## #...Citation
##   Guangchuang Yu, Li-Gen Wang, Yanyan Han and Qing-Yu He.
##   clusterProfiler: an R package for comparing biological themes among
##   gene clusters. OMICS: A Journal of Integrative Biology 2012,
##   16(5):284-287

goMFevening <- enrichGO(gene = resSig$evening_exercise$entrez,
                        universe = res$run_evening_vs_moring[, entrez], 
                        OrgDb = org.Mm.eg.db,
                        ont = "MF",
                        readable = T)

dotplot(goMFevening)

ggsave("figures/GO MF evening enrichment.pdf", height = 297/2, width = 210, units = "mm")

No sedentary evening effects. Visualizing only the exercise evening

Dexamethasone, iso & CL target gene enrichment

Generating signatures

The source scripts are not included for brevity, the sets were generated using the GEO2R function from GEO database with the groups clarified below.

# source("dex.R")
# source("iso.R")
# source("cl.R")

load("data/CL.treat.mice.Rdata")
load("data/iso.treat.3T3.Rdata")
load("data/dex.treat.3T3.Rdata")

3T3-L1:

Iso experiment: https://www-ncbi-nlm-nih-gov.ezproxy.u-pec.fr/geo/query/acc.cgi?acc=GSE43658

Dexamethasone: https://www-ncbi-nlm-nih-gov.ezproxy.u-pec.fr/geo/query/acc.cgi?acc=GSE62635 (selected only 24 hour treatment and control)

Mouse: beta3 CL experiment: https://www-ncbi-nlm-nih-gov.ezproxy.u-pec.fr/geo/query/acc.cgi?acc=GSE55934

Analyzing sets

design <- model.matrix(~0 + super_group, data = dataGroups)
colnames(design) <- gsub(colnames(design), pattern = "super_group", replacement = "")

pharmacStim <- list(iso = iso3T3[iso3T3$logFC>0 & iso3T3$adj.P.Val<0.05,"external_gene_name"], 
                    CL = ClMice[ClMice$logFC>0 & ClMice$adj.P.Val<0.05,"external_gene_name"], 
                    dex = dex3T3[dex3T3$logFC>0 & dex3T3$adj.P.Val<0.05, "Gene.symbol"])
pharmacStim <- lapply(pharmacStim, \(x){
  which(dataDGE$genes$Symbol %in% x)
})

resPharmacStim <- lapply(colnames(contrasts), function(x){
  y <- fry(y = counts, index = pharmacStim, design = design, contrast = contrasts[,x]) #|> as.data.table()
  y$contrast <- x
  y$set <- rownames(y)
  as.data.table(y)
})
resPharmacStim <- Reduce(rbind, resPharmacStim)

resPharmacStimDir <- lapply(colnames(contrasts), \(i){
  qlf <- glmQLFTest(fit, contrast=contrasts[,i])
  lapply(pharmacStim, \(ii){
    topTags(qlf[ii,], n = Inf, sort.by = "none")$table[,"logFC"]
  })
})
names(resPharmacStimDir) <- colnames(contrasts)

Plotting

dataPlot <- reshape2::melt(resPharmacStimDir) |> as.data.table()
colnames(dataPlot) <- c("logFC", "set", "contrast")

setkey(resPharmacStim, set, contrast)
dataPlot <- cbind(resPharmacStim[.(dataPlot$set, dataPlot$contrast)], dataPlot[,"logFC"])
dataPlot$FDRlog <- -log10(dataPlot$FDR)
dataPlot[FDR>0.05, FDRlog := NA]
dataPlot[FDR>0.05, FDR := NA]

dataPlot$set <- as.factor(dataPlot$set)
levels(dataPlot$set) <- c("CL316", "Dex", "Iso")
dataPlot$set <- relevel(dataPlot$set, "Dex")

dataPlot$contrast <- as.factor(dataPlot$contrast)
levels(dataPlot$contrast) <- c("Active phase\nexercise",
                               "Rest phase\nexercise",
                               "Time of day\neffect exercise",
                               "Time of day\neffect sedentary")

ggplot(dataPlot, aes(x = set, y = logFC, fill = FDR)) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_boxplot(outlier.size = 0.5) +
  scale_fill_distiller() +
  facet_wrap(~contrast, strip.position = "bottom", ncol = 4) +
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.placement = "outside",
        strip.text = element_text(size = 12),
        axis.text = element_text(size = 10, color = "black")) +
  xlab("")

ggsave("figures/Pharmacological enrichment.png", height = 297/3, width = 210, units = "mm")

Session info

sessionInfo()

## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] beepr_1.3             clipr_0.8.0           ReactomePA_1.36.0    
##  [4] clusterProfiler_4.0.5 UpSetR_1.4.0          reshape2_1.4.4       
##  [7] ggsci_2.9             patchwork_1.1.1       ggplot2_3.3.5        
## [10] openxlsx_4.2.5        org.Mm.eg.db_3.13.0   org.Hs.eg.db_3.13.0  
## [13] AnnotationDbi_1.54.1  IRanges_2.26.0        S4Vectors_0.30.2     
## [16] Biobase_2.52.0        BiocGenerics_0.38.0   dplyr_1.0.8          
## [19] data.table_1.14.2     edgeR_3.34.1          limma_3.48.3         
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.4.1        shadowtext_0.1.1       fastmatch_1.1-3       
##   [4] systemfonts_1.0.4      plyr_1.8.6             igraph_1.2.11         
##   [7] lazyeval_0.2.2         splines_4.1.0          crosstalk_1.2.0       
##  [10] BiocParallel_1.26.2    usethis_2.1.5          GenomeInfoDb_1.28.4   
##  [13] digest_0.6.29          yulab.utils_0.0.4      htmltools_0.5.2.9000  
##  [16] GOSemSim_2.18.1        viridis_0.6.2          GO.db_3.13.0          
##  [19] fansi_1.0.2            checkmate_2.0.0        magrittr_2.0.2        
##  [22] memoise_2.0.1          remotes_2.4.2          Biostrings_2.60.2     
##  [25] graphlayouts_0.8.0     enrichplot_1.12.3      prettyunits_1.1.1     
##  [28] colorspace_2.0-3       rappdirs_0.3.3         blob_1.2.2            
##  [31] ggrepel_0.9.1          textshaping_0.3.6      xfun_0.30             
##  [34] callr_3.7.0            crayon_1.5.0           RCurl_1.98-1.6        
##  [37] jsonlite_1.8.0         graph_1.70.0           scatterpie_0.1.7      
##  [40] ape_5.6-2              glue_1.6.2             polyclip_1.10-0       
##  [43] gtable_0.3.0           zlibbioc_1.38.0        XVector_0.32.0        
##  [46] graphite_1.38.0        pkgbuild_1.3.1         scales_1.1.1          
##  [49] DOSE_3.18.3            DBI_1.1.2              Rcpp_1.0.8.3          
##  [52] viridisLite_0.4.0      gridGraphics_0.5-1     tidytree_0.3.9        
##  [55] reactome.db_1.76.0     bit_4.0.4              DT_0.21               
##  [58] htmlwidgets_1.5.4      httr_1.4.2             fgsea_1.18.0          
##  [61] RColorBrewer_1.1-2     ellipsis_0.3.2         pkgconfig_2.0.3       
##  [64] farver_2.1.0           sass_0.4.0             locfit_1.5-9.5        
##  [67] utf8_1.2.2             labeling_0.4.2         ggplotify_0.1.0       
##  [70] tidyselect_1.1.2       rlang_1.0.2            munsell_0.5.0         
##  [73] tools_4.1.0            cachem_1.0.6           downloader_0.4        
##  [76] cli_3.2.0              audio_0.1-10           generics_0.1.2        
##  [79] RSQLite_2.2.10         devtools_2.4.3         evaluate_0.15         
##  [82] stringr_1.4.0          fastmap_1.1.0          ragg_1.2.2            
##  [85] yaml_2.3.5             ggtree_3.0.4           processx_3.5.2        
##  [88] knitr_1.37             bit64_4.0.5            fs_1.5.2              
##  [91] tidygraph_1.2.0        zip_2.2.0              purrr_0.3.4           
##  [94] easypackages_0.1.0     KEGGREST_1.32.0        ggraph_2.0.5          
##  [97] nlme_3.1-155           aplot_0.1.2            DO.db_2.9             
## [100] brio_1.1.3             compiler_4.1.0         rstudioapi_0.13       
## [103] beeswarm_0.4.0         png_0.1-7              testthat_3.1.2        
## [106] treeio_1.16.2          statmod_1.4.36         tibble_3.1.6          
## [109] tweenr_1.0.2           bslib_0.3.1            stringi_1.7.6         
## [112] highr_0.9              ps_1.6.0               desc_1.4.1            
## [115] lattice_0.20-45        Matrix_1.4-0           vctrs_0.3.8           
## [118] pillar_1.7.0           lifecycle_1.0.1        jquerylib_0.1.4       
## [121] cowplot_1.1.1          bitops_1.0-7           qvalue_2.24.0         
## [124] R6_2.5.1               gridExtra_2.3          vipor_0.4.5           
## [127] sessioninfo_1.2.2      MASS_7.3-55            assertthat_0.2.1      
## [130] pkgload_1.2.4          rprojroot_2.0.2        withr_2.5.0           
## [133] GenomeInfoDbData_1.2.6 grid_4.1.0             ggfun_0.0.5           
## [136] tidyr_1.2.0            rmarkdown_2.13         ggforce_0.3.3         
## [139] ggbeeswarm_0.6.0